The statistical mechanics of networks 
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We study the family of network models derived by requiring the expected properties of a graph 
ensemble to match a given set of measurements of a real-world network, while maximizing the 
entropy of the ensemble. Models of this type play the same role in the study of networks as is played 
by the Boltzmann distribution in classical statistical mechanics; they offer the best prediction of 
network properties subject to the constraints imposed by a given set of observations. We give exact 
solutions of models within this class that incorporate arbitrary degree distributions and arbitrary 
but independent edge probabilities. We also discuss some more complex examples with correlated 
edges that can be solved approximately or exactly by adapting various familiar methods, including 
mean-field theory, perturbation theory, and saddle-point expansions. 



o 

CD 



CO 

03 



i 

C 

O 
o 



> 

\o 
in 
in 
o 
^i- 
o 



-a 
c 

o 
o 



X 



I. INTRODUCTION 

The last few years have seen the publication of a large 
volume of work in the physics literature on networks 
of various kinds, particularly computer and information 
networks like the Internet and world wide web, biologi- 
cal networks such as food webs and metabolic networks, 
and social networks 0, 0, S 0| ■ This work has been di- 
vided between empirical studies of the structure of partic- 
ular networks and theoretical studies focused largely on 
the creation of mathematical and computational models. 
The construction of network models is the topic of this 
paper. 

Models of networks can help us to understand the im- 
portant features of network structure and the interplay 
of structure with processes that take place on networks, 
such as the flow of traffic on the Internet or the spread 
of a disease over a social network. Most network mod- 
els studied in the physics community arc of a practical 
sort. Typically one wishes to create a network that dis- 
plays some feature or features observed in empirical stud- 
ies. The principal approach is to list possible mechanisms 
that might be responsible for creating those features and 
then make a model incorporating some or all of those 
mechanisms. One then either examines the networks pro- 
duced by the model for rewarding similarity to the real- 
world systems they are supposed to mimic, or uses them 
as a substrate for further modeling, for example of traffic 
flow or disease spread. Classic examples of models of this 
kind are the small-world model |5j and the many differ- 
ent preferential attachment models 0,0,@j which model 
network transitivity and power-law degree distributions 
respectively. 

However, there is another possible approach to the 
modeling of networks, which has been pursued compar- 
atively little so far. An instructive analogy can be made 
here with theories of gases. There are (at least) two dif- 
ferent general theories of the properties of gases. Kinetic 
theory explicitly models collections of individual atoms, 
their motions and collisions, and attempts to calculate 
overall properties of the resulting system from basic me- 
chanical principles. Pressure, for instance, is calculated 



from the mean momentum transfered to the walls of a 
container by bombarding atoms. Kinetic theory is well 
motivated, easy to understand, and makes good sense 
to physicists and laymen alike. However, kinetic theory 
rapidly becomes complex and difficult to use if we at- 
tempt to make it realistic by the inclusion of accurate 
intermolecular potentials and similar features. In prac- 
tice, kinetic theory models either make only rather rough 
and uncontrolled predictions, or they rely on large-scale 
computer simulation to achieve accuracy. 

If one wants a good calculational tool for studying the 
properties of gases, therefore, one does not use kinetic 
theory. Instead, one uses statistical mechanics. Although 
certainly less intuitive, statistical mechanics is based on 
rigorous probabilistic arguments and gives accurate and 
reliable answers for an enormous range of problems, in- 
cluding many, such as problems concerning solids, for 
which kinetic theory is inapplicable. Equilibrium statisti- 
cal mechanics provides a general framework for reasoning 
and a powerful calculational tool for very many problems 
in statistical physics. 

Here we argue that the current commonly used mod- 
els of networks are akin to kinetic theory. They posit 
plausible mechanisms or dynamics, and produce results 
in qualitative agreement with reality, at least in some 
respects. They are easy to understand and give us good 
physical insight. However, like kinetic theory, they do not 
make quantitatively accurate predictions and provide no 
overall framework for modeling, each network model in- 
stead concentrating on explaining one or a few features 
of the system of interest. 

In this paper we discuss exponential random graphs, 
which are to networks as statistical mechanics is to the 
study of gases — a well-founded general theory with true 
predictive power. These advantages come at a price: ex- 
ponential random graphs are both mathematically and 
conceptually sophisticated, and their understanding de- 
mands some effort of the reader. We believe this effort 
to be more than worthwhile, however. Theoretical tech- 
niques based on solid statistical foundations and capable 
of quantitative predictions have been of extraordinary 
value in the study of fluid, solid state, and other physical 
systems, and there is no reason to think they will be any 
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less valuable for networks. 

We are by no means the first authors to study ex- 
ponential random graphs, although our approach is dif- 
ferent from that taken by others. Exponential random 
graphs were first proposed in the early 1980s by Hol- 
land and Lcinhardt l9j, building on statistical founda- 
tions laid by Besag |ldj |. Substantial further develop- 
ments were made by Frank and Strauss 0, 0, fl3| . 
and continued to be made by others throughout the 
1990s n~E| . In recent years a number of physicists, 
including ourselves, have made theoretical studies of spe- 
cific cases [H E3 El El IH |2l|. Today, exponential 
random graphs are in common use within the statistics 
and social network analysis communities as a practical 
tool for modeling networks and several standard com- 
puter tools are available for simulating and manipulating 
them, including Prepstar, ERGM, and Siena |22l |. 

In this paper we aim to do a number of things. First, 
we place exponential random graph models on a firm 
physical foundation, showing that they can be derived 
from first principles using maximum entropy arguments. 
In doing so, we argue that these models are not merely 
an ad hoc formulation studied primarily for their mathe- 
matical convenience, but a true and correct extension of 
the statistical mechanics of Boltzmann and Gibbs to the 
network world. 

Second, we take an almost entirely analytic approach 
in our work, by contrast with the numerical simulations 
that form the core of most previous studies. We show 
that the analytic techniques of equilibrium statistical me- 
chanics are ideally suited to the study of these models 
and can shed much light on their structure and behav- 
ior. Throughout the paper we give numerous examples 
of specific models that are solvable either exactly or ap- 
proximately, including several that have a long history 
in network analysis. Nonetheless, the particular exam- 
ples studied in this paper form only a tiny fraction of the 
possibilities offered by this class of models. There are 
many intriguing avenues for future research on exponen- 
tial random graphs that are open for exploration, and we 
highlight a number of these throughout the paper. 



II. EXPONENTIAL RANDOM GRAPHS 

The typical scenario addressed in the creation of a net- 
work model is this: one has measurements of a number 
of network properties for a real-world network or net- 
works, such as number of vertices or edges, vertex de- 
grees, clustering coefficients, correlation functions, and 
so forth, and one wishes to make a model network that 
has the same or similar values of these properties. For 
instance, one might find that a network has a degree se- 
quence with a power-law distribution and wish to create 
a model network that shows the same power law. Or one 
might measure a high clustering coefficient in a network 
and wish to build a model network with similarly high 
clustering. 



Essentially all models considered in modern work, and 
indeed as far back as the 1950s and 1960s, have been 
ensemble models, meaning that a model is defined to be 
not a single network, but a probability distribution over 
many possible networks. We adopt this approach here as 
well. Our goal will be to choose a probability distribution 
such that networks that are a better fit to observed char- 
acteristics are accorded higher probability in the model. 

Consider a set of graphs. One can use any set Sf, but 
in most of the work described in this paper J# will be the 
set of all simple graphs without self-loops on n vertices. 
(A simple graph is a graph having at most a single edge 
between any pair of vertices. A self- loop is an edge that 
connects a vertex to itself.) Certainly there are many 
other possible choices and we consider some of the others 
briefly in Sections IIII Dl and IIII Fl The graphs can also 
be either directed or undirected and we consider both in 
this paper, although most of our time will be spent on 
the undirected case. 

Suppose we have a collection of graph observables {xi}, 
i = 1 ... 7', that we have measured in empirical observa- 
tion of some real-world network or networks of interest 
to us. We will, for the sake of generality, assume that 
we have an estimate (xi) of the expectation value of each 
observable. In practice it is often the case that we have 
only one measurement of an observable. For instance, 
we have only one Internet, and hence only one measure- 
ment of the clustering coefficient of the Internet. In that 
case, however, our best estimate of the expectation value 
of the clustering coefficient is simply equal to the one 
measurement that we have. 

Let G 6 be a graph in our set of graphs and let P{G) 
be the probability of that graph within our ensemble. We 
would like to choose P(G) so that the expectation value 
of each of our graph observables {xi} within that distri- 
bution is equal to its observed value, but this is a vastly 
underdetermined problem in most cases; the number of 
degrees of freedom in the definition of the probability dis- 
tribution is huge compared to the number of constraints 
imposed by our observations. Problems of this type how- 
ever are commonplace in statistical physics and we know 
well how to deal with them. The best choice of probabil- 
ity distribution, in a sense that we will make precise in a 
moment, is the one that maximizes the Gibbs entropy 

S = -£p(G)]nP(G), (1) 

subject to the constraints 

J2P(G)x i (G) = (x i ), (2) 

G 

plus the normalization condition 

£P(G) = 1. (3) 

G 

Here Xi(G) is the value of Xi in graph G. 
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Introducing Lagrange multipliers a, {9i}, we then find 
that the maximum entropy is achieved for the distribu- 
tion satisfying 

+ afl-J]P(G) N 

G 



dP{G) 



+ X>((s i )-X;P(G)z<(G ! ) 



G 



for all graphs G. This gives 

In P(G) + 1 + a + ^2 &i x i i G ) = °> 



(4) 



(5) 



or equivalently 



P(G) 



-H{G) 



where H(G) is the graph Hamiltonian 
H(G) = ^2e iXi (G) 



(6) 



(7) 



(8) 



and Z is the partition function 
Z = c a+1 = £\ 

G 

Equations © to JSJl define the exponential random graph 
model. The exponential random graph is the distribution 
over a specified set of graphs that maximizes the entropy 
subject to the known constraints. It is also the exact 
analogue for graphs of the Boltzmann distribution of a 
physical system over its microstates at finite tempera- 
ture. 

Using the exponential random graph model involves 
performing averages over the probability distribution JBJ). 
The expected value of any graph property x within the 
model is simply 



c)=J2P(G)x(G). 



(9) 



The exponential random graph, like all such maximum 
entropy ensembles, gives the best prediction of an un- 
known quantity x, given a set of known quantities, 
Eq. 0. In this precise sense, the exponential random 
graph is the best ensemble model we can construct for a 
network given a particular set of observations. 

In many cases we may not need to perform the sum (0 ; 
often we need only perform the partition function sum, 
Eq. (JSJ), and the values of other sums can then be de- 
duced by taking appropriate derivatives. Just as in con- 
ventional equilibrium statistical mechanics, however, per- 
forming even the partition function sum analytically may 
not be easy. Indeed in some cases it may not be possible 
at all, in which case one may have to turn to Monte Carlo 
simulation, to which the model lends itself admirably. As 
we show in this paper however, there are a variety of tools 
one can employ to get exact or approximate analytic so- 
lutions in cases of interest, including mean- field theory, 
algebraic transformations, and diagrammatic perturba- 
tion theory. 



III. SIMPLE EXAMPLES 

Before delving into the more complicated calculations, 
let us illustrate the use of exponential random graphs 
with some simple examples. 



A. Random graphs 

Consider first what is perhaps the simplest of exponen- 
tial random graphs, at least for the case of fixed number 
of vertices n considered here. 

Suppose we know only the expect number of edges (m) 
that our network should have. In that case the Hamilto- 
nian takes the simple form 



H(G) = 0m(G). 



(10) 



We can think of the parameter 9 as either a field coupling 
to the number of edges, or alternatively as an inverse 
temperature. 

Let us evaluate the partition function for this Hamil- 
tonian for the case of an ensemble of simple undirected 
graphs on n vertices without self-loops. We define the 
adjacency matrix <j to be the symmetric n x n matrix 
with elements 



1 if i is connected to j, 
otherwise. 



(11) 



Then the number of edges is m 
partition function is 



Si<j a ij' and tne 



G 



i<j fif=0 i<j 



(12) 



It is convenient to define the free energy 

F = -laZ, (13) 

which in this case is 



F l • e "). 



(14) 



(Note that the free energy is extensive not in the num- 
ber of vertices n, but in the number (™) of pairs of ver- 
tices, since this is the number of degrees of freedom in 
the model.) Then, for instance, the expected number of 
edges in the model is 



G 

n\ 1 



I dZ OF 

z~de ~ ~m 



2 j e e + 1 



(15) 
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Conventionally we re-express the parameter 9 in terms of Then the partition function is 



P = 



c e + l 



(16) 



so that (m) = Q)p- 

The probability P(G) of a graph in this ensemble can 
be written 



P(G) 



[1 



I (5) 



= p m (l-p)(S)-" 1 . (17) 



In other words, P{G) is simply the probability for a graph 
in which each of the Q) possible edges appears with in- 
dependent probability p. 

This model is known as the Bernoulli random graph, 
or often just the random graph, and was introduced, 
in a completely different fashion, by Solomonoff and 
Rapoport [23( in 1951 and later famously studied by 
Erdos and Renyi 0, E3- Today it is one of the best 
studied of graph models, although, as many authors have 
pointed out, it is not a good model of most real- world net- 
works 0, One way in which its inadequacy shows, 
and one that has been emphasized heavily in networks 
research in the last few years, is its degree distribution. 
Since each edge in the model appears with independent 
probability p, the degree of a vertex, i.e., the number of 
edges attached to that vertex, follows a binomial distri- 
bution, or a Poisson distribution in the limit of large n. 
Most real-world networks however have degree distribu- 
tions that are far from Poissonian, typically being highly 
right-skewed, with a small proportion of vertices hav- 
ing very high degree. Some of the most interesting net- 
works, including the Internet and the world wide web, 
appear to have degree distributions that follow a power 
law 0, |2(|, • In the next section we discuss what hap- 
pens when we incorporate observations like these into our 
models. 



B. Generalized random graphs 

Suppose then that rather than just measuring the total 
number of edges in a network, we measure the degrees of 
all the vertices. Let us denote by ki the degree of ver- 
tex i. The complete set {ki} is called the degree sequence 
of the network. Note that we do not need to specify in- 
dependently the number of edges m in the network, since 
m = ^ Yli ki f° r an undirected graph. 

The exponential random graph model appropriate to 
this set of observations is the model having Hamiltonian 



(18) 



where we now have one parameter 9i for each vertex i. 
Noting that ki = Y]j (Xij , this can also be written 



(19) 



Z = £ expf-^T^ 
and the free energy is 



i<j &ij=0 

(20) 



F=-^ln(l+e-^+^). (21) 

i<j 

More generally we could specify a Hamiltonian 

H = '£e ij a ij , (22) 

i<j 

with a separate parameter 0^ coupling to each edge • 
Then 

z=n(i4 

i<3 



»**), F = -X>(l + e- e *>). (23) 

i<j 



This allows us for example to calculate the probability of 
occurrence pij of an edge between vertices i and j: 



Pij = (cy) 



dF 



e e 'j + 1 



(24) 



The model of Eq. I|18(l is the special case in which 
©ij = 0i+0j and the normal (Bernoulli) random graph of 
Eq. (|12|) corresponds to the case in which the parameters 
@ij are all equal. 

Sometimes it is convenient to specify not a degree se- 
quence but a probability distribution over vertex degrees. 
This can be achieved by specifying an equivalent distribu- 
tion over the parameters 0i in (|18fl . Let us define p{&) d9 
to be the probability that the parameter 9 for a vertex 
lies in the range to 6 + d9. Then, averaging over the 
disorder so introduced, the free energy, Eq. H21|) , becomes 



F = - f p{0 1 )d9 l . . . p{9 n ) d9 n ]T ln(l + e 



ln(l + c- (e+e '^p(9)p(9') d9 d9'. (25) 



The part of this free energy due to a single vertex with 
field parameter 9 is 



1 SF 



= -(n - 1) / ln(l + c- {e+e ">)p{9') d9', (26) 



n 8p{9) 



and the expected degree of vertex i with field 9i is the 
derivative of this with respect to 9, evaluated at 9f. 



(k^ = -(n-1) 
= (n-1) 



±Jln(l + c-^)p(9>)d9> 
p(9') d9' 



e e i+ e< + i ■ 



(27) 
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By a judicious choice of p(0) we can then produce the 
desired degree distribution. (See also Sec. 1111 El ) Wc 
studied this model in a previous paper 0, as a model 
for degree correlations in the Internet and other networks. 

We could alternatively specify a probability distribu- 
tion p(9) for the parameters 8y in (|22J) that couple to 
individual edges. Or, taking the developments a step fur- 
ther, one could define joint distributions for the O^ on 
different edges, thereby introducing correlations of quite 
general kinds between the edges in the model. There are 
enormous possibilities to be explored in this regard, but 
we pass over them for now, our interests in the present 
paper lying in other directions. 

One can calculate many other properties of our models. 
For example, for the model of Eq. . one can calculate 
the expectation value of any product of vertex degrees 
from an appropriate derivative of the partition function: 



(kikj 



d d 

Wide] 



z. 



(28) 



Such derivatives are correlation functions of degrees 
within the model. Similarly, derivatives of the free en- 
ergy give the connected correlation functions: 



dF 

d 2 F 
dOidOj 
d 3 F 

de~de~dOi 



(h) , (29a) 
(k i k j ) c = (k i k j )-(k i )(k j ), (29b) 

(kikjki) c 

= (hkjki) - (kikj) c (hi) - (kjki) c (h) 
- (k l k l ) c (kj) - (ki) (kj) (h) , (29c) 

and so forth. 

For instance, the two-vertex connected correlation is 



(kikj) 




for i ^ j, 
for i = j. 



(30) 



For the case of the Bernoulli random graph, which has all 
9i equal, this gives (kikj) = p(l — p) for i ^ j, where we 
have made use of Eq. ijToj) . Thus the degrees of vertices 
in the random graph are in general positively correlated. 
One can understand this as an effect of the one edge 
that potentially connects the two vertices i and j. The 
presence or absence of this edge introduces a correlation 
between the two degrees. (For a sparse graph, in which 
p = 0(n _1 ), the correlation disappears in the limit of 
large graph size.) 

In order to measure some quantities within exponen- 
tial random graph models, it may be necessary to in- 
troduce additional terms into the Hamiltonian. For in- 
stance, to find the expectation value of the clustering 
coefficient C [5j, one would like to evaluate 



(C) 



J2g C(G)e 
Z 



-H 



(31) 



which we can do by introducing an extra term linear in 
the clustering coefficient in the Hamiltonian. To measure 
clustering in the network of Eq. Ijl8jl . for example, we 
could define 



Then 



dF 



(32) 



(33) 



7=0 



Thus it is important, even in the simplest of cases, to be 
able to solve more general models, and much of the rest 
of the paper is devoted to the development of techniques 
to do this. 



C. Directed graphs 

Before we look at more complicated Hamiltonians, let 
us look briefly at what happens if we change the graph 
set over which our sums are performed. The first case 
we examine is that of directed graphs. We define to 
be the set of all simple loopless directed graphs, which is 
parameterized by the asymmetric adjacency matrix 



if there is an edge from j to i, , 
otherwise. ^ ' 



Thus, for instance, the Hamiltonian H = 9m gives rise 
to a partition function 



t^j <?ij=0 



[1 + e-*] 



2(2) 



(35) 



and a corresponding free energy. 

The directed equivalent of the more general model of 
Eq. H18J) in which we can control the degree of each ver- 
tex is a model that now has two separate parameters 
for each vertex, Of 1 and 0f ut , that couple to the in- and 
out-degrees: 

# = E(^c+c t fcr t )- ( 36 ) 

i 

Then the partition function and free energy are 

z = n( i+c ~ (c+en ) ( 3? ) 

F = -Y,H l + Q ~ [ef+er) )- (38) 

From these we can calculate the expected in- and out- 
degree of a vertex: 



1 



Mi 

dF 
89° ut 



= E 

3(¥=i) 



,(0f+0? ut ) 



i ' + 1 



e (<?f +e? ut ) + i ' 



(39) 
(40) 
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We note that J^i (^ n ) = (^° ut ), as must be the case 
for all directed graphs, since every edge on such a graph 
must both start and end at exactly one vertex. 

We can also define a probability distribution 
p(8 m ,8 ont ) for the fields on the vertices, and the devel- 
opments generalize Eqs. I|25H27I) in a natural fashion. 

We give a more complex example of a directed graph 
model in Section IIV C II where we derive a solution to 
the reciprocity model of Holland and Leinhardt using 
perturbative methods. 



D. Fermionic and bosonic graphs 

It will by now have occurred to many readers that re- 
sults like Eqs. i|21|) and (|27() bear a similarity to cor- 
responding results from traditional statistical mechanics 
for systems of non- interacting fermions. We can look 
upon the edges in our networks as being like particles in 
a quantum gas and pairs of vertices as being like single- 
particle states. Simple graphs then correspond to the 
case in which each single-particle state can be occupied 
by at most one particle, so it should come as no surprise 
that the results look similar to a system obeying the Pauli 
exclusion principle. 

Not all networks need have only a single edge between 
any pair of vertices. Some can have multiple edges or 
multiedges. The world wide web is an example — there 
can be and frequently is more than one link from one page 
to another. The Internet, airline networks, metabolic 
networks, neural networks, citation networks, and col- 
laboration networks are other examples of networks that 
can exhibit multiedges. There is no problem generalizing 
our exponential random graphs to this case and, as we 
might expect, it gives rise to a formalism that resembles 
the theory of bosons. 

Let us define our set of graphs & to be the set of all 
undirected graphs with any number of edges between any 
pair of vertices (but still no self-edges, although there is 
no reason in principle why these cannot be included as 
well). Taking for example the Hamiltonian, Eq. (|22J) . and 
generalizing the adjacency matrix, Eq. so that <Jij 

is now equal to the number of edges between i and j, we 
have 



n-ij between vertices i and j , which is given by 



{<?ij} i<j i<3 o-» 3 '=0 



KJ 



and 



(41) 



(42) 



i<3 



The equivalent of the probability p^ of an edge appearing 
in the fermionic case is now the expected number of edges 



n ij — 



1 



dF 



<)(-),, e e - 1 ' 



(43) 



Note that this quantity diverges if we allow Gjj — * 0, 
a phenomenon related to Bose-Einstein condensation in 
ordinary Bose gases. 

For the special cases of Eqs. ()18|l and (|1C)|> . we have 



^ = £ln(l-e-^+^), 

i<j 



*« ~ e e,+e 3 _ ! 



and 



F = 



^ ~~ n0 



e° -1 



(44) 



(45) 



respectively. The connected correlation between the de- 
grees of any two vertices in the latter case is 



{kikj) 



d 2 F 



do 2 (e e -i) 2 ' 



(46) 



for i =/= j. Thus the degrees are again positively correlated 
and the correlation diverges as 9 — > 0. 



E. The sparse or classical limit 

In most real-world networks the number of edges m is 
quite small. Typically m is of the same order as n, rather 
than being of order n 2 . Such graphs are said to be sparse. 
(One possible exception is food webs, which appear to be 
dense, having m = 0(n 2 ) H^-) The probability of an 
edge appearing between any particular vertex pair 
is of order 1/n in such networks. Thus, for example, in 
the fermionic case of the network described by the Hamil- 
tonian J22J, Eq. (|2~l|l tells us that e Qij must be of order n 
in a sparse graph. The same is also true for the bosonic 
network of the previous section. This allows us to ap- 
proximate many of our expressions by ignoring terms of 
order 1 by comparison with terms of order e ij . We refer 
to such approximations as the "sparse limit" or the "clas- 
sical limit," the latter by analogy with the corresponding 
phenomenon in quantum gases at low density. 

In particular, the equivalent of Eq. (|24|l for either 
fermionic or bosonic graphs in the classical limit is p^ — 



For the case of Eq. i|18[l. it is 



Pij 



(47) 



so that each edge appears with a probability that is a sim- 
ple product of "fugacities" e~ 6i defined on each vertex. 
The classical limit of this model has been studied previ- 
ously by a number of other authors [13, 123, kill i 
although again developed and justified in a different way 
from our presentation here; generally the edge probabil- 
ity l|47|l has been taken as an assumption, rather than a 
derived result. 



7 



For a given distribution p(9) of 9, the expected degree 
of a vertex, Eq. itTTT . is 

{ki) = (n-l)e- 8i Je- e 'p(e')d0', (48) 

which is simply proportional to e~ 6i . So we can produce 
any desired degree distribution by choosing the corre- 
sponding distribution for 9. 

F. Fixed edge counts 

Another possible choice of graph set is the set of 
graphs with both a fixed number of vertices n and a fixed 
number of edges m. Models of this kind have been ex- 
amined occasionally in the literature .16] and, if we once 
more adopt the view of the edges in a graph as particles, 
they can be considered to be the canonical ensemble of 
network models, where the variable edge-count models of 
previous sections are the grand canonical ensemble. As 
in conventional statistical mechanics, the grand ensemble 
is often simpler to work with than the canonical one, but 
progress can be made sometimes be made in the canoni- 
cal case by performing the sum over all graphs regardless 
of edge count and introducing a Kronecker <5-symbol into 
the partition function to impose the edge constraint: 

Z = ^2S(m,m(G))e- H , (49) 

G 

where m is the desired number of edges. 

For instance, the fixed edge-count version of the gen- 
eralized random graph, Eq. I|22|). would be one in which 




(50) 

where we have made use of the integral representation 
for the (^-function 

6(m,m)= [ e 27d( ™~ m),? d v . (51) 
Jo 

The sum over graphs is now in the form of the partition 
function for the grand canonical version of the model, 
but with 0jj — > Qij + 27ri?7, giving the field parameters 
an imaginary part. Thus, from Eq. (|23|) 

Z= [ d77e 2 ^TT(l + e- (e ^ +27ri, ' ) ). (52) 

In general the integral cannot be done in closed form, 
which is why fixed edge-count graphs — and canonical en- 
sembles in general — are avoided. The integral can in 



principle be carried out term by term for any finite n, 
but doing so is tantamount to performing the sum over 
all graphs with fh edges explicitly, so there is little to be 
gained by the exercise. 

It is also possible to have a bosonic graph with a fixed 
number of edges — one would simply sum over the set of 
graphs that have fh edges with any number of them being 
permitted to fall between any given pair of vertices. 

We will not discuss further either fixed edge-counts or 
bosonic networks in this paper, concentrating instead on 
the grand canonical fermionic ones, which are more useful 
overall. However, essentially all of the results reported in 
the remainder of the paper can be generalized, with a 
little work, to these other cases if necessary. 



IV. MORE COMPLEX HAMILTONIANS 

Outside of the models described in the previous sec- 
tions, and some minor variations on them, we know of 
few other exponential random graph models that are ex- 
actly solvable. (One exception is the reciprocity model of 
Holland and Leinhardt for which we derive an exact 
solution in Sec. IIV C 1\ ) To make further progress one 
must turn to approximate methods. There are (at least) 
three types of techniques that can yield approximate ana- 
lytic solutions for exponential random graph models. The 
first and simplest is mean-field theory, which works well 
in many cases because of the intrinsically high dimen- 
sionality of network models; usually these models have 
an effective dimensionality that increases with the num- 
ber of vertices n, so that the thermodynamic limit of 
n — > oo also corresponds to the high dimension limit in 
which mean-field theory becomes accurate. Nonetheless, 
there are many quantities, such as those depending on 
fluctuations, about which mean-field theory says noth- 
ing, and for these other methods are needed. In some 
cases one can use non-perturbative approaches based on 
the Hubbard-Stratonovich transform or similar integral 
transforms, which are very effective and accurate but 
suitable only for models with Hamiltonians of specific 
forms polynomial in the adjacency matrix. More gen- 
erally, one can use perturbation theory, which may in- 
volve larger approximations (although they are usually 
well controlled), but is applicable to Hamiltonians of es- 
sentially any form. 

We discuss all of these approaches here. As an example 
of their application, we use one of the oldest and best- 
studied of exponential random graphs, the 2-star model. 
The Hamiltonian for the 2-star model is 

H = 9m — as, (53) 

where m is the number of edges in the network and s is 
the number of "2-stars." A 2-star is two edges connected 
to a common vertex. (The minus sign in front of the 
parameter a is introduced for later convenience.) 

The quantities m and s can be rewritten in terms of 
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the degree sequence thus: 



2 ^ ] ki? s — 2 ^ y ki(ki 1). 



(54) 



Substituting these expressions into Eq. I|53|> . we can 
rewrite the Hamiltonian as 



(55) 



where J = \{n— l)a and B — — ^(9+a). (The factor (n— 
1) in the definition of J is also introduced for convenience 
later on.) 

Noticing once again that ki = J"^- , where the vari- 
ables CTij are the elements of the adjacency matrix, we 
can also write 



H 



J 



- <?ij <Tik -B^cxij. (56) 

ijk ij 



We study the 2-star model in the fermionic case in which 
each vertex pair can be connected by at most a single 
edge, and within the grand canonical ensemble where the 
total number of edges is not fixed. Generalization to the 
other cases described above is of course possible, if not 
always easy. 



A. Mean- field theory 

The variables <7ij can be thought of as Ising spins resid- 
ing on the edges of a fully connected graph, and hence the 
2-star model can be thought of as an Ising model on the 
edge-dual graph of the fully connected graph plj. (The 
edge-dual G* of a graph G is the graph in which each edge 
in G is replaced by a vertex in G* and two vertices in G* 
are connected by an edge if the corresponding edges in 
G share a vertex.) Using this equivalence, the mean-field 
theory of the 2-star model can be developed in exactly 
the same way as for the more familiar lattice-based Ising 
model. 

We begin by writing out all terms in Eq. I|56|) that 
involve a particular spin cry : 



H((7ij) = -CFi 



J 



1 ^ 

k 



crjk + cr k j) + 2B 



where we have explicitly taken account of all the ways 
in which o~ij can enter the first term in the Hamiltonian. 
(We have also dropped the term 2 Ju^jin — 1) required to 
correctly count the terms diagonal in , since it vanishes 
in the large n limit.) 

Then, in classic mean-field fashion, we approximate 
the local field by its average: 



J 



1 ^ 

k 



(o lk + a kl + a jk + a kj ) + 2B -> 4 Jp + 2B, (58) 



J= 1.5 



J= 1 



7 = 0.5 



0.5 § 



-1 

external field B 



FIG. 1: The mean-field solution for the connectance p — 
(k) /(n — 1) in the 2-star model from Eq. 1601 . for values of 
the coupling J below, at, and above the phase transition. For 
the case J = 1.5 we are in the symmetry broken phase and 
the hysteresis loop corresponding to the high- and low-density 
phases of the system is clearly visible. 



where, as before, p = (a) is the mean probability of an 
edge between any pair of vertices, which is also called 
the connectance of the graph. Then H(o~ij) = — (4 Jp + 
2B)aij, and we can write a self-consistency condition for 
p of the form 



P 



-H(o 



l—p g-i?(CT i3 =0) 

Rearranging, this then gives us 

c 4./p+2S 



_ e 4.Jp+2B_ 



V 



I _|_ e 4Jp+2B 



\ [tanh(2 Jp + B) + i\. 



(59) 



(60) 



For J < 1 this equation has only one solution, but 
for J > 1 there may either be one solution or, if B is 
sufficiently close to — J, there may be three, of which the 
outer two are stable. Thus when B is close to — J we have 
a bifurcation at J c = 1, a continuous phase transition to 
a symmetry broken state with two phases, one of high 
density and one of low. We show in Fig. a plot of the 
solution of l|6UII which displays clearly the characteristic 
hysteresis loop of the symmetry broken state. 

Along the "symmetric line" B = — J there is always a 
solution p = i to Eq. I|6U|) (although it may be unstable), 
and along this line we can think of p ~ \ as a standard 
order parameter for the model which is zero in the high- 
symmetry phase and non-zero in the symmetry-broken 
phase. We can define a critical exponent (3 in the usual 
fashion by 



VP- 



{J -I) 



(61) 



as J — > + , giving j3 = i, which is the usual Ising mean- 
field value and should come as no surprise, given the 
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equivalence mentioned above between the 2-star model 
and the Ising model. One can define other critical expo- 
nents as well, which are also found to take Ising mean- 
field values. For instance, as we showed in the vari- 
ance x of the connectance, which plays the role of a sus- 
ceptibility, goes as x ~ I J ~ 1|~ 7 m the vicinity of the 
phase transition with 7 = 1. 



B. Non-perturbative approaches 

We can go beyond the mean-field approximation of the 
previous section by making use of techniques borrowed 
from many-body theory. The developments of this sec- 
tion follow closely the lines of our previous paper on this 
topic |18| , and rather than duplicate material needlessly, 
the reader is referred to that paper for details of the calcu- 
lation. Here we merely summarize the important results. 

The evaluation of the partition function for the 2-star 
model involves a sum of terms of the form e . The study 
of interacting quantum systems has taught us that such 
sums can be performed using the Hubbard-Stratonovich 
transform. We start by noting the well-known result for 
the Gaussian integral: 



-atj) 2 



Making the substitutions a — ► (n — 1) J and 
ki/(n — 1), and rearranging, this becomes 



(62) 



3 Jk?/(n-l) 



(n-l)J 



X / 

-oo 



-(n-l)J0|+2J0ifci 



(63) 



Then the partition function is 
J 



^exp 

G 

(n-l)J 



n — 

n/2 



J 2)0 exp(-(n-l)j£# 



£expf£(2J<& + B)k, 



(64) 



where we have interchanged the order of sum and inte- 
gral. 

The sum over graphs now has precisely the form of the 
partition function sum for the model of Eq. 118(1. and 
from Eq. I|2U|) we can thus immediately write down the 
partition function 



-j*r(4>) 



(65) 



where the quantity 

Jt>{4>) = (n-l)J^^I^ln(l + c 2 *+«+ 2B ) 

i ijtj 

- inln((n- 1)J), (66) 



is called the effective Hamiltonian. 

Thus we have completed the partition function sum 
for the 2-star model, but at the expense of introducing 
the auxiliary fields {4>i} which must be integrated out to 
complete the calculation. The integral cannot, as far as 
we are aware, be evaluated exactly in closed form but, 
as we showed in .18], it can be evaluated approximately 
using a saddle-point expansion, with the result that the 
free energy of the 2-star model is given to leading order 
in the expansion by 



F = n(n - l)J<t>l - |n(n - 1) ln(l + < 



where 



+ |(n-l)ln(l-2J0o(l-0o)), 



6 = i[tanh(2J0 o + B) +1] 



(67) 



(68) 



is the position of the saddle point, i.e., the maximum of 
the Hamiltonian on the real-(/> line. 

Note that Eq. (|68|l is identical to the mean-field equa- 
tion, Eq. (|60l) . for the connectance p of the 2-star model. 
Thus, 0o is the connectance of the model within the 
mean-field approximation and the saddle-point expan- 
sion, as is typically the case in such calculations, is an 
expansion about the mean-field solution. 

From the free energy we can derive a number of quan- 
tities of interest. We showed in [T||, for instance, that 
the variance of vertex degree in the model is given by 



(k 2 )-(kf = (n-l) 



0o(l - <j>o) 
1- 2J^,(1 -<f>oY 



(69) 



which has a gradient discontinuity but no divergence at 
the phase transition. (This quantity is, by contrast, zero 
everywhere within mean- field theory.) 



C. Perturbation theory 

Exponential random graphs also lend themselves nat- 
urally to treatment using perturbation theory. Here we 
describe the simplest such theory, which is roughly equiv- 
alent to the high-temperature expansions of conventional 
thermal statistical mechanics. Expansions of this t ype 
have been examined previously by Burda et al. [lOL l2(l ] 
for Strauss's model of a transitive network [l2lll8j|. Here 
we develop the theory further for general exponential ran- 
dom graphs. 

The fundamental idea of perturbation theory for ran- 
dom graphs is the same as for other perturbative meth- 
ods: we expand about a solvable model in powers of the 
coupling parameters Qi in the Hamiltonian. We write the 
Hamiltonian for the full model in the form H = Hq + Hi, 
where Ho is the Hamiltonian for the solvable model and 
Hi takes whatever form is necessary to give the correct 
expression for H. Then the partition function is p^ll20| 



(70) 
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where Zq = J^g e ~ H ° is the partition function for the un- 
perturbed Hamiltonian, and (■ ■ ■ )o indicates an ensemble 
average in the unperturbed model. 

The only case that has been investigated in any de- 
tail is the one where we expand around a random graph, 
Hq = dm, so that the averages in Eq. Ij70|l are averages 
in the ensemble of the random graph. (It is possible for 
9 to be zero, so this choice for Hq does not place any re- 
striction on the form of the overall Hamiltonian. If 8 = 
then the expansion is precisely equivalent to an ordinary 
high-temperature series.) However, for Hamiltonians H 
that give significant probability to networks substantially 
different from random graphs, the perturbation theory 
cannot be expected to give accurate answers at low or- 
der. In theory there is no reason why one could not ex- 
pand about some other solvable case, although no such 
calculations have been done as far as we are aware. One 
obvious possibility, which we do not pursue here, is to ex- 
pand around one of the generalized random graph forms, 
Eqs. JED and 

Typically, to make progress with Eq. (JJDJ, we will ex- 
pand the exponential in a power series of the form 

fc=0 

In practice, Hi normally contains a coupling constant, 
such as the constant J in the 2-star model of Eq. I|55|l . and 
thus our expression for the perturbed partition function 
is an expansion in powers of the coupling. 

In this section, we apply the perturbation method to 
two example models. First, we study a simple model pro- 
posed about a quarter of a century ago by Holland and 
Leinhardt which is exactly solvable by this method. 
Then we illustrate the application of the method to the 
2-star model and compare its performance against the 
approximate saddle-point expansion results of the previ- 
ous section. 

1. The reciprocity model 

Our first example of perturbation theory is a directed 
graph model. In the real world, many directed graphs 
display the phenomenon of reciprocity: a directed edge 
running from vertex A to vertex B predisposes the net- 
work to have an edge running from B to A as well. Put 
another way, the network has a higher fraction of vertex 
pairs that are joined in both directions than one would 
expect on the basis of chance ( "mutual dyads" in the par- 
lance of social network analysis). Behavior of this kind is 
seen, for example, in the world wide web, email networks, 
and neural and metabolic networks |33l l34l l35j| . 

Holland and Leinhardt Q have proposed a exponential 
random graph model of reciprocity, which we study here 
in a simplified version. As we now show, the perturba- 
tion expansion for this model can be written down to all 
orders and resummed to give an exact expression for the 
partition function. 



The Hamiltonian for the model is 

H = H o + Hi=0m- ar, (72) 

where m is the total number of (directed) edges in the 
graph, and r is the number of vertex pairs with edges run- 
ning between them in both directions. The unperturbed 
Hamiltonian Hq is that of an undirected random graph 
f Sec. IIII H|| with partition function given by Eq. (|35[1 and 
each directed edge present with independent probability 
p = (c 6 + 1) . The perturbation Hi can be written in 
terms of the adjacency matrix (|34f> as 

Hi = -ar = -a)aij<jji. (73) 

Then the perturbation series, Eq. (J7TJ, for the full Hamil- 
tonian is 

<7 °° k 

I -£*<**>.■ < 74 » 

u fc=0 

with 

( rfc )o = zZ ■ ■ ■ Yl ■ ■ ■ °iH3kVi k ik) ■ (75) 

il <jl ik <jk 

Thus the partition function is written as an expansion in 
powers of a whose coefficients are correlation functions 
of elements of the adjacency matrix, calculated within 
the ordinary random graph. If we can evaluate these 
correlation functions, at least up to some finite order, we 
can also evaluate the perturbed partition function. 

Since all edges are present or absent independently of 
one another in the random graph, the correlation func- 
tions factor: 

(fl2Cr2l(7340'43)o = (°"12)o (°"2l)o ( 034)0 (°"43)o = P* 1 

and so forth. The only exception is in cases where two 
or more of the elements cry being averaged are the same. 
In that case, noting that a™- = cry for any n, we have 
results like 

(<7l 2 CT2ieri2CT2l) = (eri 2 ) (c2l) = p 2 ■ (77) 

To evaluate expressions such as (I75)) . therefore, we need 
to count the number of independent elements cry that ap- 
pear in each term. This can be difficult for some models, 
but for the reciprocity model it is quite straightforward. 
The question we need to answer is this: if we choose k 
pairs of vertices from the possible pairs, with 

duplication allowed, how many ways are there of choos- 
ing them such that exactly q pairs will be distinct? Each 
such way makes a contribution of p 2q to the partition 
function J73J. 

Let CLj^q be the number of ways of choosing the pairs 
such at a particular set of q distinct pairs are chosen at 
least once each. Note that ak,i = 1 for all k. Then, from 
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Eqs. (EU) and JTSJ) 



Z 



9=1 



^ )<7,(<*)p 2 * 



where the function 



9q 



(*) = £ 



fe=l 



it! 



(78) 



(79) 



is the exponential generating function for afc i9 . 

Now the number of ways of choosing k pairs such that 
all choices are made from a particular set of size q, but 
without the constraint that each pair in the set appear 
at least once, is just q k . Thus 
equivalently 



5Ll = <? fc , or 



Ojfe,, 



9-1 

E 



Multiplying by z k /kl and summing over fc = 1 . 
gives 



ff? (z) 



a? 2 



9-1 

E 



(80) 
. oo, this 

(81) 



which immediately implies that g q (z) = (e z — l) 9 , by 
induction on 1811) with the initial condition gi(z) — 

E^i * Vi A ! = e * - !■ 

Substituting this result into Eq. (|78() then gives us our 
solution: 



Z 



(2) //;;.. 

E 

q=l 



2 ; )(e a -l)V 9 = [l + (c Q -iy 



I (2) 



(82) 



Or, making use of Eqs. IT^jl and (TTC?|l 

'l + (e Q -l)p 2 1^ 



Z = 



1-p 



(83) 



and F = — In Z in the normal fashion. 

From these expressions we can, for instance, obtain the 
mean number of edges (to) and the mean number (r) of 
pairs of vertices connected by edges running both ways 
from 

dF dF dF 

<■»>-« -*»-i>a>. W— < 8 ") 

A quantity of interest in directed networks is the reci- 
procity |3J] , which is the fraction of edges that are recip- 
rocated. This quantity is found to be on the order of tens 
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FIG. 2: Fhe reciprocity and connectance of the model of Hol- 
land and Leinhardt [9( for p = 0.01. The solid lines represent 
the exact solution, Eq. I|8^. and the points are Monte Carlo 
simulation results for systems of n = 1000 vertices. 



of percent in networks such as the world wide web. The 
reciprocity for the model of Holland and Leinhardt is 



2(r) 



pe u 



(85) 



(to) 1 — p + pe a 

In Fig. |3 we show the reciprocity, along with the con- 
nectance of the network, as a function of a for the case 
p = 0.01. There is no phase transition or other unex- 
pected behavior in this model: the measured properties 
are smooth functions of the independent parameters. No- 
tice that there is a substantial range of values of a over 
which the connectance is low and the graph realistically 
sparse, but the reciprocity is still high, with values simi- 
lar to those seen in real networks. 



2. Example 2: The 2-star model 

As our second example of the application of pertur- 
bation theory, we return to the 2-star model introduced 
at the beginning of Sec. IIVI Unlike the case of the reci- 
procity model in the preceding section, perturbation the- 
ory does not lead to an exact solution of the 2-star model 
but, as we now show, we can get an approximate solution 
by studying the perturbation expansion to finite order — a 
different approximation from the saddle-point expansion 
of Sec. Ell 

We divide the Hamiltonian H = 9m— as into an unper- 
turbed part Hq = 9m, which is the normal Bernoulli ran- 
dom graph, and a perturbation Hamiltonian H\ = —as. 
Then, following Eq. I|71|) . the partition function for the 
full model is given by 



(86) 



1=0 



12 

W V <W ^LIV 



The number of 2-stars is 

i j<k 

and therefore 

iljl<fcl ii,ji<ki 

Our strategy is to evaluate the series Ij86|) up to some 
finite order in a to get an approximate solution for Z, but 
there is a problem. Each term in the series corresponds 
to states of the graph that have the corresponding num- 
ber of 2-stars: the term in (s) , for instance, counts the 
number of graphs that have a 2-star in any position in 
the graph. This is not enough for our purposes however. 
Realistic graphs will have not a finite number but a finite 
density of 2-stars in them, and the number of such graphs 
is counted by terms that appear at infinite order in the 
perturbation expansion in the limit n — > oo. So, without 
going to infinite order as we did in the reciprocity model, 
we are never going to get meaningful results from our 
expansion. 

Similar problems appear in ordinary statistical me- 
chanics and the solution is well known. Instead of ex- 
panding the partition function, we form an expansion for 
the free energy. We can write the free energy as 

F=-]nZ = -\nZ -]n—=F +Fi, (89) 

Zo 

where i*b is the free energy of the unperturbed network 
and F\ = —Ih(Z/Zq). Now we expand F\ as a power 
series in a of the form 

Fx = -ah-^h-^h-..., (90) 

where we have made use of the fact that F\ — when 
a = 0. Substituting into Z/Z = e~ Fl , we get 

Y q = l+afi + ^(f2+fi) + ^(f3+ZhA+fi)+0(a% 

(91) 

and comparing terms with Eq. Ij86(l . we find 

fx = (s) , (92a) 
h = {-s 2 ) Q - fl (92b) 
h = (A- Wi-A 3 > (92c) 

and so forth. These are the cumulants of s within the en- 
semble defined by the unperturbed network. If we expand 
s in the form of Eq. 1|87[) then they are connected correla- 
tions of elements of the adjacency matrix — "connected" 
because individual elements of the adjacency matrix are 
uncorrelated, so that all terms in the cumulants vanish 
unless they involve sets of 2-stars that share one or more 
edges. (Note that sharing a vertex, as in the more fa- 
miliar spin models of traditional statistical mechanics, is 



FIG. 3: The diagrams contributing to the first three orders 
in the perturbation expansion of the free energy of the 2-star 
model in powers of a. 

not a sufficient condition for being connected. The fun- 
damental degrees of freedom in a network are the edges.) 

We will proceed then as follows. We calculate the free 
energy Fi in terms of connected correlations up to some 
finite order in a and from this we calculate the partition 
function Z = Z^e~ Fl . Even though F\ is known only 
to finite order, our expression for Z will include terms 
with all powers of the connected correlations in it, via 
the expansion of the exponential, and hence will include 
graphs with not only a finite number but a finite den- 
sity of 2-stars. This idea, which will be routine for those 
familiar with conventional diagrammatic many-body the- 
ory, is entirely general and can be applied to any model, 
not just the 2-star model. In essence, the series given by 
c~ Fl is a partial resummation to all orders of the parti- 
tion function, including some but not all of the contribu- 
tions to Z from disconnected correlations of arbitrarily 
high order. 

Let us see how the calculation proceeds for the case 
of the 2-star model, to order a 3 , as above. The leading 
0(a) term in Fi is simple: 

/1 = ( s )o = E E °*>o = n ( U 2 V ' (93) 

i j<k » ' 

Since we are primarily interested in large networks, we 
can approximate this expression by its value to leading 
order in n, which is in 3 p 2 . 

The second term, at order a 2 , is more complicated 
because there are several different ways in which two 2- 
stars may combine to share one or more edges. In order to 
keep track of these different contributions, we make use of 
a diagrammatic representation similar to that emp loyed 
by Burda et al. for Strauss's transitivity model |19|. Fig- 
ure|3ji shows the single diagram contributing to fx , which 
gives the result in Eq. I|93|l . Figure Ob shows the three 
diagrams that contribute to ji- It is an assumption of 
our notation that each edge that appears in a diagram 
is distinct. Thus the third diagram in Fig. |3|d, which 
represents the case in which the two 2-stars fall on top 
of one another, must be depicted separately, rather than 
being considered a special case of the first diagram. This 
turns out to be a good idea, since this term has a differ- 
ent functional form from the first diagram, and neither 
diagram is necessarily negligible by comparison with the 
other. 
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In general the basic "Feynman rules" for interpreting 
the diagrams are: 

1. each edge contributes a factor of p; 

2. each vertex contributes a factor of n; 

3. the numerical multiplier is the number of distinct 
ways in which the diagram can be decomposed into 
overlapping 2-stars such that each edge occurs at 
least once, divided by the symmetry factor for the 
diagram. (The symmetry factor is the number of 
distinct permutations of the vertices that leave the 
diagram unchanged.) 

Then for the connected correlation functions one must 
subtract all other ways of composing lower order dia- 
grams to make the given diagram, as in Eq. I|92|) . 

To see how these rules work in practice, let us apply 
them to the first diagram in Fig. ^jp. This diagram has 
four vertices and three edges, which gives a factor of n 4 p 3 , 
by the first two rules. The diagram can be decomposed 
into two 2-stars in 6 different ways, but the symmetry 
factor is also 6, so we end up with n 4 p 3 x 6/6 = n 4 p 3 . 
The contribution to the diagram from the term — ff in 
Eq. (|92bp is —n 4 p 4 , so the final value of the diagram 
is n 4 (p 3 — p 4 ) to leading order in n. Proceeding in a 
similar fashion, the other diagrams of Fig. |3Jd contribute 
n 4 (p 3 — p 4 ) and ^n 3 (p 2 — p 4 ), respectively. The dia- 
grams for the 0(a 3 ) term are shown in Fig.^, and are 
more complicated, but routine to evaluate using the rules 
above. The final expressions for the fs are: 
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FIG. 4: The connectance of the 2-star model calculated from 
the saddle-point expansion of Sec. llVfjl fsolid line), and from 
the first- (dotted line), second- (dashed line), and third- 
order (dot-dashed line) perturbation expansions. The calcu- 
lations were performed along the symmetric line B — — J of 
Sec. II V Al where the half-filled graph with connectance ^ is 
always a solution of the mean-field equation 16811 . For J > 1 
there exist two symmetry-equivalent stable solutions in addi- 
tion to the half-filled graph. We show only the sparser of the 
two. Inset: the density of 2-stars in the same model. 



W P 2 , 



h 

h = \n 3 {l-p)p 2 {l + Anp), 

h 



(94a) 
(94b) 



§n 3 (l - p)p 2 (l + Unp + 29n 2 p 2 - 43n 2 p 3 )- 



(94c) 



Note that we have retained the leading order terms in n 
separately at each order in p, since we have no knowledge 
a priori about the relative magnitude of n and p. In a 
sparse graph, we expect that p will be of order 1/n, in 
which case it may be possible to neglect some terms. 

Once we have the expansion of F±, it is straightfor- 
ward to calculate statistical averages from derivatives of 
the free energy in the normal fashion. For example, the 
expected number of 2-stars in the network is given by 



dF 
da 



8Fi 
da 



{s)=-— = -^ = h + af 2 + \a z h + 0(a 3 ). (95) 
And the expected number of edges is 



{m) =m =p{p - l \-dp- + ^p- 

= \n 2 p + n 3 (l -p)p 2 a[l + 5(1 + 6np - 8np 2 )a 



+ i(l + 21njJ + 58nV 



2„3 



180nV 



129n W 2 ] . 

(96) 



In Fig. 21 we show the connectance 2 (to) /n 2 and the 
density of 2-stars 2 (s) / n 3 calculated from the saddle- 
point method of Sec. IIVBI and from the expressions 
above, at first, second, and third order. As the fig- 
ure shows, the perturbation expansion agrees with the 
non-perturbative method at high and low values of J = 
i(n — l)a, and markedly better for the third-order ap- 
proximation than for the first- and second-order ones. 
However, in the region of the phase transition at J c = 1 
the agreement is poor, as we would expect. In this region 
there will be large critical fluctuations and hence contri- 
butions to the free energy from large connected diagrams 
that are entirely missing from our series expansion. Pre- 
sumably by extending the perturbation series we can de- 
rive successively more accurate answers in the critical re- 
gion. We also note that the perturbation expansion gives 
results only for the sparse phase in the symmetry-broken 
region. 

We have here studied in detail two examples of the 
treatment of exponential random graphs by perturba- 
tion theory (and another can be found in Ref. 19]). The 
techniques we have used, however, are entirely general 
and diagrammatic theories similar to these, with simi- 
larly simple "Feynman rules," can be derived for other 
examples as well. 
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V. CONCLUSIONS 

In this paper we have discussed exponential random 
graphs, which in both a figurative and a quantitative 
sense play the role of a Boltzmann ensemble for the study 
of networks. Exponential random graphs are a formally 
well-founded framework for making predictions about the 
expected properties of networks given specific measure- 
ments of properties of those networks. We have shown in 
this paper how they can be derived in moderately rigor- 
ous fashion from maximum entropy assumptions about 
probability distributions over graph ensembles. 

We have given many examples of particular calcula- 
tions using exponential random graphs, starting with 
simple random graph models that have linear Hamil- 
tonians, many of which have been presented previously 
by other authors, albeit it with rather different motiva- 
tion. In most cases these linear models can be solved 
exactly, meaning that we can derive the partition func- 
tion or equivalently the free energy of the graph ensemble 
exactly in the limit of large system size. 

For nonlinear Hamiltonians it appears possible to find 
exact solutions only rarely, but we have been able to find 
approximation solutions in several cases using a number 
of different methods. Taking the particular example of 
the 2-star model, we have shown how its behavior can be 
understood using mean-field theory, perturbation theory, 



and non-perturbative methods based on the Hubbard- 
Stratonovich transform. We have also given one example, 
the reciprocity model of Holland and Leinhardt, that is 
exactly solvable by evaluating its perturbation expansion 
to all orders. 

The results presented in this paper are only a tiny 
fraction of what can be done with exponential random 
graphs. There are many interesting challenges, both 
practical and mathematical, posed by this class of mod- 
els. Exploration of the behavior and predictions of spe- 
cific models as functions of their free parameters, devel- 
opment of other approximate solution methods, or ex- 
pansion of those presented here, and the development of 
models to study network phenomena of particular inter- 
est, such as vertex-vertex correlations, effects of hidden 
variables, effects of degree distributions, and transitivity, 
are all excellent directions for further research. We hope 
to see some of these topics pursued in the near future. 
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